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Abstract. We propose and analyze a finite element method for a semi- 
stationary Stokes system modeling compressible fluid flow subject to a Navier- 
slip boundary condition. The velocity (momentum) equation is approximated 
by a mixed finite element method using the lowest order Nedelec spaces of the 
first kind. The continuity equation is approximated by a standard piecewise 
constant upwind discontinuous Galcrkin scheme. Our main result states that 
the numerical method converges to a weak solution. The convergence proof 
consists of two main steps: (i) To establish strong spatial compactness of the 
velocity field, which is intricate since the element spa<;es are only div or curl 
conforming, (ii) To prove that the discontinuous Galerkin approximations 
converge strongly, which is required in view of the nonlinear pressure func- 
tion. Tools involved in the analysis include a higher integrability estimate for 
the discontinuous Galerkin approximations, a discrete equation for the effec- 
tive viscous fiux, and various renormalized formulations of the discontinuous 
Galerkin scheme. 
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1. Introduction 

The purpose of this paper is to prove convergence of a finite element method for 
the semi-stationary barotropic compressible Stokes system 

Qt + div{gu) = 0, in (0,T) x Q, (1.1) 

-fiAu- XDdiyu + Dp{q) = f, in(0,T)xf2, (1.2) 

with initial data 

Q\t=o = Qo, on CI. (1.3) 

Here O is a simply connected, bounded, open, polygonal domain in {N = 2,3), 
with Lipschitz boundary dfl, and T > is a fixed final time. The unknowns are 
the density g = Q{t, a;) > and the velocity u = u{t, x) S K^, with x G Cl and 
t e (0, T). We denote by div and D the usual spatial divergence and gradient 
operators and by A the spatial Laplace operator. 

The pressure function is assumed to be of the form p{g) = ag^, with a > 
(Boyle's law). Typical values of 7 ranges from a maximum of | for monoatomic 
gases, through | for diatomic gases including air, to lower values close to 1 for 
polyatomic gases at high temperatures. Throughout this paper we will always 
assume that 7 > 1. The case 7=1 can also be treated; indeed, it is simpler 
since the pressure function is linear. Furthermore, the viscosity coefiicients jjL, A are 
assumed to be constant and satisfy /i > 0, A^A + 2^ > 0. 

The study of the system (1.1)-(1.2) can be motivated in several ways. Firstly, 
the system can be used as a model equation for the barotropic compressible Navier- 
Stokes equations. This might be a reasonable approximation for strongly viscous 
fiuids, where convection may be neglected. Secondly, Lions [16] use solutions of 
(1.1)-(1.2) to construct solutions to the barotropic compressible Navier-Stokes 
equations. 

Among many others, the semi-stationary system (1.1)-(1.3) has been studied 
by Lions in [16, Section 8.2]. He proves the existence of weak solutions and some 
higher regularity results. In particular, weak solutions was proven to be unique in 
the case of periodic boundary conditions or when the equations are solved on the 
hole of R^. Uniqueness was not obtained in the case of regular Dirichlet boundary 
conditions and moreover higher regularity results was only shown to hold locally. 

In this paper we impose the following boundary conditions, which axe relevant 
in the context of geophysical fiuids and shallow water models: 

u-y = Q, on{0,T)xdn, (1.4) 

and 

curlu = 0, on (0,T) X an if JV = 2, 

(15) 

cmluxu = 0, on (0,T) X an if iV = 3, 

where u denotes the unit outward normal to dfl. The first condition is a natural 
condition of impermeability type on the normal velocity. The second condition is in 
the literature often referred to as the Navier-slip condition. It can be interpreted as 
a viscous dissipation term at the boundary (more precisely "non-dissipation" since 
this term is equal to zero) [16]. 

In some geophysical applications, conditions like (1.4)-(1.5) are preferred over 
the classical Dirichlet condition since the latter necessitates expensive calculations 
of boundary layers. Of more importance to this paper, the boundary conditions 
(1.4)-(1.5) will allow us to use the finite element method in a solution space that 
can be split into two orthogonal parts in terms of a discrete version of the Hodge 
decomposition, a fact that will play a crucial role in our analysis. 
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Although many numerical methods have been proposed for the compressible 
Stokes and Navier-Stokes equations, the convergence properties of these methods 
are mostly unsettled, especially in several spatial dimensions. Ultimately, it is not 
clear if these numerical methods converge to a weak solution as the discretization 
parameters tend to zero. In one dimension, the available results are due to Hoff and 
his collaborators [23, 24, 25]. All these results apply to the compressible Navier- 
Stokes equations in Lagrangian coordinates, and moreover require the initial density 
to be of bounded variation. Interesting results regarding the existence and long time 
behavior of solutions to the one dimensional compressible Navier- Stokes have also 
been obtained using semi-discrete finite difference schemes in [13, 14, 4], again in 
Lagrangian coordinates with the initial density of bounded total variation. In more 
than one spatial dimension, we refer to a recent paper [11] in which a convergent 
numerical method for a stationary compressible Stokes system is proposed. The 
Stokes system considered in [11] is similar to (1.1)-(1.2) with linear pressure and 
no temporal dependence. 

Let us now discuss our choice of numerical method for the semi-stationary Stokes 
system. For the discretization of (1.1) we utilize a discontinuous Galerkin scheme 
based on piecewise constant approximations in space and time. The discontinuous 
Galerkin scheme was introduced more than 30 years ago [15, 20] and has since then 
undergone a blooming development, cf. [5, 6] for a review. In the context of linear 
transport equations with rough (i.e., non-Lipschitz) coefficients, a discontinuous 
Galerkin scheme, with piecewise polynomial approximations of arbitrary degree in 
the spatial variable and piecewise constant or linear approximations in the temporal 
variable, has recently been analyzed by Walkington in [22] . The work [22] is further 
developed in [17] for the variable-density incompressible Navier-Stokes equations. 

Let us now turn to the velocity (or momentum) equation (1.2). By introducing 
the vorticity w = curlu as an auxiliary unknown, keeping in mind the vector 
identity —A = curl curl— Ddiv, we can recast the momentum equation as 

^curlw — (A + /i)£'divM + Dp{g) = f, (L6) 

where we suppress the time variable t (we refer the reader to subsequent sections 
for more precision). Hence the velocity equation (1.2), together with the boundary 
conditions (1.4) (1.5), admits a formulation that lends itself naturally to a mixed 
finite element method [12, 18, 19]. 

Denote by WjJ*'^'^ the vector fields u on fl for which div u G L'^ and u-uIqq = 0, 
and by Wg"""^'^ the vector fields w onQ for which curl w £ and w x = 0. We 
choose corresponding mixed finite element spaces Vh C Wq^^''^ and Wh C Wq^'^^''^ 
based on Nedelec's elements of the first kind [18]. The mixed finite element method 
seeks functions Wh € W/, and Uh&Vh such that 

/ jjL curl WhVh + [(a* + •^) div Uh - p{Qh)] div Vh dx = / fhVh dx, 
Jn Jn 

WhVh - cmlrihUh dx = 0, 



n 

for all {r]h, Vh) S Wh x Vh, where Qh, fh are given piecewise constant functions. 

Let us denote the numerical solution of the semi-stationary Stokes system by 
{Qh,Wh,Uh) = {Qh,Wh,Uh){t,x). The main goal is to prove that {{gh,Wh,Uh)}h>o 
converges to a weak solution, at least along a subsequence. The challenging issue 
is to show that the density approximations Qh, which on the outset is only weakly 
compact in L^, in fact converges strongly. Strong convergence is mandatory if we 
want to recover the semi-stationary Stokes system when taking the limit in the 
discrete equations as /i — > 0. Related to this issue, the above mixed method enjoys 
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some advantages over the traditional finite element method based on elements. 
In particular, the approximation spaces Wh and Vh satisfy 

Vh = cmlWh + Zh, 

for some C Vh satisfying Zh J- curl Wh ■ An immediate consequence of this 
discrete Hodge decomposition is that upon writing Uh = curlr/^, + z^, we see that 
only Zh is coupled to the density Qh and moreover that curlw^, and hence Wh, 
only depends on the data /. More importantly, equipped with the discrete Hodge 
decomposition, we can separate the quantity Pes{Qh,Uh) = P{Qh) — (A + /x) dvvuh 
from the vorticity. The quanity Peff(£'/i) w/i) is the so-called effective viscous flux 
[16] associated with our discrete equations. The fact that we can separate the 
effective viscous flux from the vorticity makes it possible to prove the following 
weak continuity property: 

lim jj Pcs_{Qh,Uh) Qh dxdt — jj Pcff g dxdt (Poff, are weak limits), (1.7) 

which is the decisive ingredient in the proof of strong convergence of the density 
approximations Qh- Related to (1.7), we prove a higher integrability estimate on 
the pressure ensuring that p{Qh), and thus also PcsiOhiUh), is weakly compact in 
L^. The energy estimate only provides a uniform bound on p{gh) in L°°{L^), so a 
priori it is not even clear that p{gh) converges weakly to an integrable function. Our 
strong convergence argument is inspired by the work of Lions on the compressible 
Navier-Stokes equations, cf. [16]. 

As part of the analysis, we also show that QhUh converges weakly to gu, where g 
and u arc weak limits of gh and Uh, respectively. This convergence is not immediate 
since the element spaces utilized for the velocity approximations are merely div or 
curl conforming. In view of the discrete continuity equation (discontinuous Galerkin 
scheme), we easily obtain a bound on {gh)i in, say, 

Li(pi/-i.i). To conclude we 

need a spatial translation estimate of the form 

Ilit/i - M;j(-, • + 0IIl2(£,2) ^ as 1^1 0, uniformly in /i. (1.8) 

In view of the discrete Hodge decomposition, we will actually only need (1.8) for 
weakly curl free approximations with a L"^ bounded divergence. 

For velocity fields that are independent of time t, (1.8) implies the compact- 
ness of {uh}h>o- time independent case, it is known that weakly curl free 
approximations with L"^ bounded divergence is compact in provided the approx- 
imation spaces satisfy the commuting diagram property [7]. However, despite the 
fact that the element spaces used here satisfy this property, the inclusion of time 
in Uh{t, x) makes earlier results inadequate. Specifically, to apply known result we 
would need L°° control in time of the velocity approximations. Unfortunately, this 
is not available in general for our problem. As a consequence, we shall provide a 
direct argument for the spatial translation estimate (1.8). 

We wish to point out that although the boundary conditions (1.4) (1.5) are 
not covered by Lions' results [16], his proofs can be adapted to yield existence, 
uniqueness, and regularity results for (1.1)-(1.2) with the boundary condtions (1.4)- 
(1.5). We will not pursue this project here, except for the existence part, which will 
be an immediate consequence of our convergence result. However, let us remark 
that the Navier-slip condition (1.5) is technically easier to handle than a Dirichlet 
condition, both from a mathematical and numerical point of view. The primary 
reason for this lies in the need for solutions of the auxiliary problem 

divt; = /, curlv = 0. (1.9) 

If J^fdx = 0, the function v will satisfy the boundary conditions (1.4)-(1.5). In 
other situations, like periodic boundary conditions or when the equations are solved 
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on M , the boundary values of v does not matter. However, it is evident that v 
cannot be required both to satisfy Dirichlet boundary conditions and (1.9). Thus, 
(1.9) can only be required to hold locally whenever Dirichlet boundary conditions 
are imposed. To avoid "localizing" various discrete arguments, which sometimes can 
require elaborate work, we have chosen to consider the Navier-slip type condition 
(1.5) instead of the no-slip Dirichlet condition. 

This paper is organized as follows: In Section 2, we introduce notation and list 
some basic results needed for the later analysis. Moreover, we recall the usual 
notion of weak solution and introduce a mixed weak formulation of the velocity 
equation. Finally, we introduce the finite element spaces and review some of their 
basic properties. In Section 3, we present the numerical method and state our main 
convergence result. The existence of a solution to the discrete equations is confirmed 
in Section 4. Section 5 is devoted to deriving basic estimates. In Section 6, we prove 
the main convergence result stated in Section 3. The proof is divided into several 
steps (subsections), including convergence of the continuity scheme, weak continuity 
of the discrete viscous flux, strong convergence of the density approximations, and 
convergence of the velocity scheme. 

2. Preliminary material 

2.1. Some functional spaces and analysis results. We make frequent use of 

the divergence and curl operators and denote these by div and curl, respectively. In 
the 2D case, we will denote both the rotation operator taking vectors into scalars 
and the curl operator taking scalars into vectors by curl. This confusing but rather 
standard notation greatly simplifies all subsequent arguments allowing identical 
treatment of the 2D and 3D cases. 
We will also make use of the spaces 

t^div,2(^) = {t; e L\n) : divw e L'^ifl)} , 
t^'^'^'-i.2(fi) = {ve L^{Q) : cmlv e L^ifl)} , 

where v denotes the unit outward pointing normal vector on dfl. If u e W'^^'^''^{Cl) 
satisfies v ■ = 0, we write v G Wo'^'^(f2). Similarly, v e Wo"'^''^(f2) means 
V e W^^^''^{Q) and v x u\gQ = 0. In two dimensions, is a scalar function and 
the space Wo^'^''^(ri) is to be understood as Wq'^{^1). To define weak solutions, we 
shall use the space 

W{Q) = {ve i^(f2) : dvfv e L^{Q),cmlv e L'^{Q),v ■ u\eQ = 0} , 

which coincides with W;^''' '"^ (n) n W'^'^'^in). The space W{n) is equipped with 
the norm = + l|div t;||^2(s^) + !|curl7;||^2(o)- K is known that \\-\\^ is 

equivalent to the norm on the space {v G H^{Q) : v ■ vlon = O}, see, e.g., [16]. 
The space W(n) admits a unique orthogonal Hodge decomposition 

W{fl) =cmlS{fl) + DA-'^Ll{fl), (2.1) 

where S{Ct) = {v e Vr^'2(n) : curlv e W'''''^{Ct)}, is the inverse Neumann 
Laplace operator, and Lq denotes the functions on fl that have zero mean. 

For the convenience of the reader we list some basic functional analysis results 
to be used in the subsequent arguments (for proofs, see, e.g., [9]). Throughout the 
paper we use overbars to denote weak limits, with the underlying spaces being 
(silently) given by the context. 

Lemma 2.1. Let O be a bounded open subset of M.^^ , M > 1. Suppose g:M.^ 
(—00,00] is a lower semicontinuous convex function and {t'n}„>i is a sequence of 
functions on O for which Vn ^ v in L^{0), g{vn) € L^{0) for each n, g{vn) — ^ 
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g{v) in L^{0). Then g{v) < g{v) a.e. on O, g{v) G L^(0), and Jpgiv) dy < 
liminfn^oo Jo di'^n) dy. If, in addition, g is strictly convex on an open interval 
{a,b) C K and g{v) = g{v) a.e. on O, then, passing to a subsequence if necessary, 
Vn{y) v{y) for a.e. y e {y e O I v{y) G (a, 6)}. 

Let X be a Banach space and denote by X* its dual. The space X* equipped 
with the weak-* topology is denoted by X^^^^^^, while X equipped with the weak 
topology is denoted by Xwoak- By the Banach- Alaoglu theorem, a bounded ball in 
X* is (7(X*, X)-compact. If X separable, then the weak-* topology is metrizable 
on bounded sets in X*, and thus one can consider the metric space C ([0, T]; X^^^-^ 
of functions w : [0, T] — !■ X* that are continuous with respect to the weak topology. 
We have Vn ^ v m. C {[Q,T\] X^^^^) if {vnit) , 4>) x\x {v{t) , 4>) x* ,x uniformly 
with respect to t, for any (p G X. The following lemma is a consequence of the 
Arzela-Ascoli theorem: 

Lemma 2.2. Let X be a separable Banach space, and suppose w„: [0,T] X* , 
n = 1,2,..., is a sequence for which H^'nllioo^o t]-x*) — some constant C 

independent of n. Suppose the sequence [0, T] 3 t ^ {vn{t),^)x* ,x, n = 1, 2, . . . , 
is equi- continuous for every <I> that belongs to a dense subset of X . Then w„ belongs 
to C {[0,T]; X^^^-^) for every n, and there exists a function v e C {^,T]; X^^^-^ 
such that along a subsequence as n —> oo there holds — > w m C ([0, T]; X^^^^. 

In what follows, we will often obtain a priori estimates for a sequence i 

that wc write as "w„ Gb X^^ for some functional space X. What this really means 
is that wc have a bound on that is independent of n. 

2.2. Topological degree in finite dimensions. Our numerical method consti- 
tutes a nonlinear-implicit discrete problem. We will prove the existence of a solution 

to this problem by a topological degree argument [8]. 

Denote by d{F,fl,y) the Z-valued (Brouwer) degree of a continuous function 
F : Q ^ at a point y G M.^\F{dS) relative to an open and bounded set 
n CR'^. For notational convenience, let us reformulate the definition of degree 
so that it applies directly in our finite element setting. Indeed, below we define 
dsh{F,Sh,qh) with F : Sh ^ Sh being a continuous finite element mapping, Sh 
being a bounded subset of a finite element space Sh, and qh being a function in Sh- 

Definition 2.3. Let Sh be a finite element space, || • || be a norm on this space, 
and introduce the bounded set 

Sh = {qheSh;\\qh\\<C}, 

where C > is a constant. Let {(7,}^^ be a basis such that span{cri}^^ = Sh and 
define the operator Us : Sh — > by 

M 

^Bqh = {qi,q2,---,qM), qh = J^gic^z- 

i=l 

The degree ds^ {F, Sh, qh) of a continuous mapping F : Sh ^ Sh sd, qh G Sh\F{dSh) 
relative to Sh is defined as 

dsdF,Sh,qh) = d (llsF{U^'),UBSh,IlBqh) ■ 

The next lemma is a consequence of the properties of the degree d{F, fi, y), cf. [8]. 

Lemma 2.4. Fix a finite element space Sh, and let dsi^{F, Sh,, qh) be the associated 
degree of Definition 2. 3. The following properties hold: 

(1) dsh{F, Sh, qh) does not depend on the choice of basis for Sh- 

(2) ds^{ld,Sh,qh) = l. 
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(3) dsh{H{- ,a), Sh, 9/1 (a)) is independent of a € J := [0,1] for H:Sh x J ^ Sh 
continuous, Qh ■ J ^ Sh continuous, and qh{c() ^ H{dSh, a) Va e [0, 1]. 

(4) ds,{F,Sh,qh)^0 ^ F-\qh)^<D. 

2.3. Wecik and renormalized solutions. 

Definition 2.5 (Weak solutions). We say tliat a pair {g, u) of functions constitutes 
a weak solution of tiie semi-stationary compressible Stokes system (1.1) (1-2) witti 
initial data (1.3) and Navier-slip type boundary conditions (1.4)-(1.5) provided the 
following conditions hold: 

(1) {g,u) € L°°{Q,T\L'i{n)) x L^{Q,T-W{n))- 

(2) Qt + diw{Qu) = in the weak sense, i.e, G C~([0, T) x Q), 

/ / Q{(l)t + uD(l)) dxdt+ [ go(l)\t=o dx = 0; (2.2) 

(3) -iiAu-\Dd\vu + Dp{Q) = / in the weak sense, i.e, V<^ G C°°((0,T) xH) 
for which </» • i/ = on (0, T) x dQ., 

I I /ticurlucurl^ + [(/X + A) divu — p(£i)] div dxdt = / / fcpdxdt. 
Jo Jq Jo Jq 

For the convergence analysis we shall also need the DiPerna-Lions concept of 
renormalized solutions of the continuity equation. 

Definition 2.6 (Renormalized solutions). Given u E i^(0, T; yV(f2)), we say that 

g G L°°(0,T; L''{il)) is a renormalized solution of (1.1) provided 

B{g)t + div {B{g)u) + b{g) div w = in the weak sense on [0, T) x Q, 

for any B G C[0, oo) n C^{Q, oo) with B{Q) = and h{g) := gB'{g) - B{g). 

We shall need the following lemma. 

Lemma 2.7. Suppose {g,u) is a weak solution according to Definition 2.5. If 
g G L^((0,r) X f2)), then g is a renorm,alized solution according to Definition 2.6. 

Proof. Let {g, u) be a weak solution. Then u G i^(0, T; iJ^(fi)). As the boundary 
of ri is Lipschitz, the velocity field u{t) can be extended to the full space such 
that u{t)\a = u{t) and 



ll^lll,2(0,T;ffi(R")) — ||M||j;,2(o,T;If 



where u{t) denotes the extension of u{t). If we extend g{t) to by setting 
g[t) = g{t)lQ, we get 

gt + div(^ w) =0 in the weak sense on [0, T) x R^. 

Now, to conclude the proof, we appeal to a well-known lemma from [16] stating 
that the square-integrable weak solution g is also a renormalized solution. □ 

2.4. A mixed formulation. In view of the Navier-slip boundary condition (1.5), 

it is natural to introduce the vorticity w = curl it as an independent variable, 
thereby turning the velocity equation into (1.6). This immediately leads to the 
following mixed formulation, which acts as a motivation for our choice of numerical 
method: Determine functions 

{w,u) G L2(0,T; Wo""''''(l))) X i2(0,T; Wo'*'^''(f2)) 
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such that 

i-T 



I I II cmX wv +[{iJ.+ X)AiY u — p{Q)\d\YV dxdt = I I fvdxdt, 
Jo Jo, Jo Jn ^2 3^ 



WT] — curlrju dxdt = 0, 



n 

for all (r/,-y) S L^^q^ T; Wo"''^(r!)) x ^^(o. T; w^jf "'"(r!)). 

In order to arrive at the weak formulation (2.3), we have utilized the integration 
by parts formula 



/ T} curl u dx = / u curl rj dx + 
Jn Jn Jd 



u{r] X ly) dS{x). (2.4) 

an 

It follows as an immediate consequence of the Stokes Theorem and will be applied 
multiple times throughout the paper. 

The upcoming goal is to prove that a sequence of approximate solutions, denoted 
by {{gh,Wh,Uh.)}f^^Q, converge to a limit {g,w,u) satisfying (2.2) and (2.3); the 
term "converge" is made precise in a forthcoming section. Having constructed such 
a limit, it follows immediately that the pair {g,u) is a weak solution according to 
Definition 2.5, thereby completing the analysis. 

2.5. Finite element spaces and some basic results. Upon inspection of the 
spatial spaces entering the weak formulations stated above, we see that they can 
be related through a De Rham sequence. In two dimensions this reads 

W^-^ Wo'""'" Ll . 0, 

while in three dimensions the corresponding sequence is 

W^''^ Wt'^ Ll > 0. 

These sequences are exact in the sense that the null space of one operator exactly 
matches the image of the next. This perspective on the spaces is actually useful as 
we see that this is precisely how the quantities u, and g relate to each other. 

It follows from (2.3) that the vorticity w is decoupled from the density g, which 
is an important consequence of our choice boundary condition and this fact is of 
relevance to the convergence analysis. Moreover, the subsequent analysis relies 
heavily on the solvability of the problem (or more precisely a discrete version of it) 

div V = q, V ■ u = on dfl, 

for some given right-hand side qinL"^. In particular, it is important for us to extract 
from this problem some control on curlw. From the above Dc Rham sequence, we 
see immediately that there exists solution v which is weakly curl free. In the 
continuous setting this is enough to conclude that curlu = 0; indeed, the Hodge 
decomposition (2.1) combined with the fact that v is weakly curl free implies v = Ds 
for some scalar s. 

Motivated by these remarks, we shall in the next section present a numerical 
method that utilizes finite element spaces satisfying a discrete version of the above 
De Rham sequence. More precisely, we will replace Wq™'''^ and W;^"^^ by the 
lowest order Nedelec finite element spaces of the first kind (but other spaces are 
possible) with vanishing degrees of freedom at the boundary dCl. Let us denote 
these spaces by Wh and Vh respectively. It is well known that the spaces Wh,Vh 
together with the space Qh of piecewise constants (cf. the ensuing section for missing 
details) satisfies in three dimensions the following exact discrete De Rham sequence: 

Sh Wh Vh QHnLUn) > 0, 
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where Sh, is the usual scalar linear Lagrange element. In the two dimensional case 
this sequence still holds, but now the spaces Sh and Wh are equal and thus the 
sequence does not contain the gradient operator. All finite element spaces are 
defined with respect to a given tetrahedral mesh E^, of Q. 
We introduce the canonical interpolation operators: 

: n w^2,2 ^ . ^curi,2 ^ ^2,2 ^ 

using the available degrees of freedom of the involved spaces. That is, the operators 
are defined 

(nf s) (xi) = s(xi), G Afh] 

j^{ll'^w)xv dS{x) = jwxv dS{x), Veef/j; 

j (n^w) • u dS{x) = J v-f dS{x), vr e r,,; 

/ n^g dx= q dx, WE e Eh, 
Je Je 

where Th, £h, and J\fh, denote the set of faces, edges, and vertices, respectively, of 
Eh- Then it is well known that the following diagram commutes: 

<l <l <l 

Wh Vn Qh. 

Remark 2.8. The interpolation operators Ilf , , and H]^ , are defined on function 
spaces with enough regularity to ensure that the corresponding degrees of freedom 
are functional on these spaces. This is reflected in writing Wg"'^''^ fl W^'^ instead 
of merely VF"^"""''^ and so on. 

In view of the above commuting diagram, we can deflne the spaces orthogonal 
to the range of the previous operator, i.e., 

W°-^ := {wh G Wh; cml Wh = 0}^ n Wh, 

:= {Vh G Vh-, divvh = 0}^ n Vh, 

to obtain decompositions 

Wh = DSh + Wl^'^, 

Vh = cmlWh + V°'^, (2.5) 
and the discrete Poincare inequalities 

Wvhh^^n) < C||divt;,||^.(^) , yv G (2.6) 

W^hWma) < C||curl«;^||^.(^) , Vt,; G W^^'^ . (2.7) 

Thus, with this configuration of elements we are able to perform unique Hodge type 
decompositions of the discrete vector fields. As an example, we immediately have 
the existence of a function Vh G VJ^*'^ satisfying 

divVhlE = qh\E, yE G Eh, 

for any given Qh&Qh^ {/q Qh dx = O}. 

The following lemma summarizes well-known error estimates satisfied by the in- 
terpolation operators. The estimates are derived from the Bramble-Hilbert lemma 
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using scaling arguments. We however note that care must be taken when mapping 
functions in Wh and V/j to a reference element (cf. [3, 18]). 

Lemma 2.9. There exists a constant C > 0, depending only on the shape regularity 
of Eh and the size offl, such that for any 1 < p < oo, 



h - ^ML.m + h ||div(t; - n^t;) ||^,(^) < Ch^ \\D^v\\^,^^^ , r = 1, 2, 
\\w - UfwW^^^^^ + h ||curl(«; - UYw)\\^^^^^ < Ch'\\D'w\\LP(^n), s = 1,2, 

for all (f) e W^-P{n),v e W-''P{n), and w G W^'P{n). 

In what follows, we will need the following lemma. It follows from scaling argu- 
ments and the equivalence of finite dimensional norms. 

Lemma 2.10. There exists a constant C > 0, depending only on the shape regu- 
larity of Eh, such that for 1 < q,p < oo, and r = 0,1, 

for any E G Eh and all polynomial functions ^h S Fk{E), fc = 0, 1, — 

The next result follows from scaling arguments and the trace theorem. 

Lemma 2.11. Fix any E G Eh and let (j) G W^''^{E) be arbitrary. There exists a 
constant C > 0, depending only on the shape regularity of Eh such that, 

Uh^r) < Ch--- {UWl^e) + hmWi^.^E)) , yreThn dE. 

3. Numerical method and main result 

In this section we define the numerical method and the state the convergence 
theorem. The proof of this theorem is deferred to subsequent sections. 

Given a time step At > 0, wc discretize the time interval [0,T] in terms of the 
points = mAt, m — 0, . . . , AI, where we assume that MAt = T. Regarding the 
spatial discretization, we let {Eh}h be a shape regular family of tetrahedral meshes 
of n, where h is the maximal diameter. It will be a standing assumption that h 
and At are related such that At = ch, for some constant c. By shape regular we 
mean that there exists a constant k > such that every E € Eh contains a ball of 
radius Xe > where He is the diameter of E. Furthermore, we let F/j denote 
the set of faces in Eh- Throughout the paper, we will use the three dimensional 
terminology (tetrahedron, face, etc.) to denote both the three dimensional and the 
two dimensional case (triangle, edge, etc). 

On each element E G Eh, we denote by Q{E) the constants on E. The functions 
that are piecewise constant with respect to the elements of a mesh Eh are denoted 
by Qh = Qh{^)- Next, on each E G Eh, we denote by W{E) the lowest order 
space of curl-conforming Nedelec polynomials of first kind [18]. In two dimensions, 
W{E) is the space of linear scalar polynomials on E and is totally determined by 
it's value at the vertices of E. In three dimensions, each member of W{E) is of the 
form 

/..\ 

a-l- 6 X 2/ , a,6 e M^, 

and is totally determined by the following degrees of freedom: J^w ■ Te dS(x) for 
all edges (not faces) e of the element E, where Tg is the unit tangential vector on e. 
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On each element E e E^, we denote by V{E) the lowest order space of div- 
conforming Nedelec polynomials of first kind [18]. In two dimensions, it is the 
Raviart-Thomas polynomial space on E. Each member of V{E) is of the form 



a e M^5e 



and is totally determined by the following degrees of freedom: J^v ■ u dS{x), for 
all faces T of the element E, where is a unit normal vector on T. 

The element spaces Wh = W/i(ri) and Vh = V^(n) are formed on the entire 
mesh Efi by matching the degrees of freedom of the polynomial space W{E) and 
V{E), respectively, on each face F G F^. In addition, wc incorporate the boundary 
conditions by letting the degrees of freedom of the spaces Wh and Vh vanish at the 
faces on the boundary. 

Before defining our numerical method, we shall need to introduce some additional 
notation related to the discontinuous Galerkin scheme. Concerning the boundary 
dE of an element E, wc write /+ for the trace of the function / achieved from within 
the clement E and /_ for the trace of / achieved from outside E. Concerning a 
face F that is shared between two elements E_ and E^, we will write for the 
trace of / achieved from within E^ and /_ for the trace of / achieved from within 
Here i?_ and E-^- are defined such that v points from to _E+, where ly is 
fixed (throughout) as one of the two possible normal components on each face F. 
We also write [/]r = /+ — f- for the jump of / across the face F, while forward 
time-differencing of / is denoted by = — J™. To denote the set of inner 

faces of Th we will use the notation F^ = {F € F/^; F ^ d^l}. 

Let us now define our numerical method for the semi-stationary Stokes system 
(1.1)-(1.2) augmented with the boundary conditions (1.4) and (1.5) (note, however, 
that in the definition below the boundary conditions are built into the finite element 
spaces and not listed explicitly). 

Definition 3.1 (Numerical method). Let {Qhi^)} h>o ^® ^ sequence in Qh{^) that 

satisfies g1 > for each fixed h > and gi" qq a.e. in and in L^{il) as h 0. 

Set fh{t, •) = /^"'(•) := ^ n^/(s, •) ds, for t € (t„_i, t„), m = 1, . . . , M. 
Now, determine functions 

{gZ',w^^,u^^) e Qhin) X Whin) X Vhin), m=l,...,M, 

such that for all (ph G Qhi^), 

[ eT<l>h dx - At [ {e^{K-,.)+ + e':^{K-i.r)[cPh]rdS{x) 

^^^^ (3.1) 

Qh~^4'h dx, 

and for all {r]h,Vh) € Wh{n) x Vh{fl), 

[ ^l curl w]!:vh + [(m + A) div < - p{q^)] div Vh dx ^ [ fJ^Vh dx, 

w^Vh - cm\r]h dx = 0, 



(3.2) 

^"hilh - Uh cm\r]h dx = 0, 

/a 

for m = 1, . . . , M . 

In (3.1), [uh ■ v)^ = max{ti/j • v,Q} and {uh ■ v)^ = min{tt/i • v,Q}, so that 
Uh- v = {uh ■ z^)"*" + {uh ■ y)~ , i.e., in the evaluation of q{u ■ u) at the face F the 
trace of g is taken in the upwind direction. 



12 



K. H. KARLSEN AND T. K. KARPER 



Remark 3.2. Using the identity 



EeE^ JdE\dn 



rer 

we can state (3.1) on the following form: 



/ g^ct^h dx + AtJ2 f ■ + ^- « • ^)~) <t>h dS{x) 

Jq JdE\dn 

9T~'<Ph dx. 



For each fixed h > 0, the numerical solution {(^?™, ■u^/^", it™)}^=o extended to 
the whole of (0, T] x f2 by setting 

{Qh,Wh,Uh){t) = iQ'^,w^,u'^), te{tm-utm], m=l,...,M. (3.4) 

In addition, we set gh{0) = q1- 

Our main result is that, passing if necessary to a subsequence, {{Qh^Wh,, Uh)}h>o 
converges to a weak solution. More precisely, there holds 

Theorem 3.3 (Convergence). Suppose f € L'^{{0,T) x and qq G i^(f^), 7 > 1. 
Let {{gh,Wh,Uh)}f^yQ be a sequence of numerical solutions constructed according 
to (3.4) and Definition 3.1. Then, passing if necessary to a subsequence as h ^ 0, 
Wh^w in L2(0,T; Wo™''''2(J7)), Uh ^ u in L'^iO,T;W^^''^{n)), qhUh Qu in 
the sense of distributions on (0,T) x Q., and Qh ^ Q «-e. in (0, T) x Q., where the 
limit {q,w,u) satisfies the mixed formulation (2.3), and consequently {g,u) is also 
a weak solution according to Definition 2.5. 

This theorem will be an immediate consequence of the results stated and proved 
in Sections 4-6. 

4. Numerical method is well defined 

In this section we show that there exists a solution to the discrete problem given 

in Definition 3.1. However, we commence by obtaining a positive lower bound for 
the density, recalling that the approximate initial density f'JK') strictly positive. 

Lemma 4.1. Fix any m = 1,...,M, and suppose G Qh{^), € V/i(fl) 

are given bounded functions. Then the solution q^{^ e Qh^) of the discontinuous 
Galerkin scheme (3.1) satisfies 



min£.™(x) > min£-"-i(a;) ( ^ „ / ^ . 

Consequently, if Q^~^{-) > 0, then gi™(-) > 0. 

Proof Let E e Eh he such that g"^]^ < g'^lj^ ^E e Eh, and insert into (3.3) the 
test function (j)h & Qh{^), defined by 



Integrating by parts then yields 



XG E, 
otherwise. 



0h 



E = -^J ie7i< ■ + s-{uh ■ '^r) dS{x) + e^-'^ 

\E\ JdE\dn 



MIXED FEM FOR A SEMI-STATIONARY STOKES SYSTEM 13 

At f 

= -Atig^dwu^) U - ^ / . ie^- e+){uH ■ v)-dS{x) + er'U. 

-c/ JdE\dn 



ldE\dn 

>-At(^jrdivoi^+er'U' 

where we have also used the relation Uh- 1' = {uh ■ v)'^ -\- {uh • is constant 

on E, and that g'^ attains its minimal value on E. Consequently, 

n"'\ - > n"^~^ I ^ 
Qh \e - Sh I 



1 + Af||div-u'^"|U=c(o) 



□ 



We now turn to the existence of solutions to our nonlinear-implicit discrete 
problem. We will apply a topological degree argument, thereby reducing the proof 

to exhibiting a solution to a linear problem. 

Lemma 4.2. For each fixed h > 0, there exists a solution 

{g';:,w]:,u';:)€Qf,{n)xWh{n)xVn{n), p;r(-)>0, m=l,...,M, 
to the nonlinear-implicit discrete problem posed in Definition 3.1. 

Proof. We argue by induction. Assume for m = 1, . . . , fc — 1 that there exists a 
solution 

to the discrete problem of Definition 3.1. Here and below we denote by Q'^{^) the 
strictly positive functions in Qh{Q). Moreover, the norm ||-|| on Sh is defined by 

The claim is that we can find a solution for m = k: 

{glwluDeSh. (4.1) 
To this end, we introduce the mapping 

H : Q+{n) X Wh{n) X Vhin) X [0, 1] ^ Qh{n) X Wh{n) x Vh{n), 

H{gh,Wh,Uh,a) = {zh{a),yh{a),Xh{a)) , 
where the triplet {zh{a),yh{oi),Xh{cx)) is defined by 

Zh{a)(j)h dx = a y2 / {Q+{uh ■ + g-{uh ■ f)') (ph dS{x) 
Jn EeEh "'S'E;\9n 

+ l^^^^^^^h dx, y^H&Qhm, 

/ Xh{a)vh dx = iJ. curl WhVh + (A + /x) div Uh div Vh dx 
Jq Jn 

^)divvhdx- fhVhdx, yvh^Vh{^), 



ot [ P{gh) 
Jn 



/ VhVh dx= WhVh - Uhcmlrjh dx, \/r]h e W/,(f2). 
Jn Jn 

Solving H{gh,Uh,'Wh,l) = is equivalent to finding a solution (4.1) to the 
nonlinear-implicit discrete problem posed in Definition 3.1. 

Let us fix an arbitrary a £ [0, 1], and consider a solution {gh{a),Uh{cx),Wh,{a)) 
belonging to Q^{^) x W'h{^) x Vh{fl) of the corresponding problem 

H{gh{a),Uh{a),Wh{a),a) = 0. 

We claim that there is a constant > 0, independent of a, such that 

\\'Wh{a)\\wonri.2^Q^ + ||Mfc(a)|lwdi-,2(n) + lb/»(Q!)|II,(a) < Cf- (4.2) 
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Indeed, repeating the arguments leading to estimate (5.15) in Section 5 we conclude 
that (4.2) holds with = C'||/|||2(o T-L'^(a)) + Il^h~^lll7(n)' where the constant C 
is the constant appearing in (5.15). Here, wc have also used that At = 0{h). Let 
Sh C Sh = Qti'^) X Vhi^) X Wh{^) be a ball of sufficiently large radius, cf (4.2). 
Then, since every solution of H{gh, Uh, Wh, a) = lies strictly inside Sh, 

O^H{dSh,a), Vae[0,l]. (4.3) 

We claim that H{-,-) is continuous on Sh x [0,1]. Let {gh,Uh,Wh) € Sh- By 
equivalence of norms on finite dimensional spaces, the functions gh,Uh,Wh are 
bounded on In view of this and Lemma 2.10, the claim follows. 

By the virtue of (4.3) and the continuity of H{-, •), we have by Lemma 2.4 that 

{H{-, a),Sh, 0) is independent of a e [0, 1]. 

The proof will be completed by proving that ds^{H{-,a ~ 0),Sh,0) 0. To see 
this, observe that the problem H{gh, Uh, Wh, 0) = is equivalent to finding a triplet 
{gh,Uh,Wh) e Sh satisfying 



/ gh(l>hdx= / gl ^(t)hdx, e Q?,(n), (4.4) 

Jn Jn 

IJ,cmlwhVh + (A + /u) divuhdiyvh dx = fhVh dx, Vvh e Vh{0.), 



j 

Jn 



(4.5) 

cmlrjhUh - WhVh dx = 0, Vjy/i e Wh{fl). 

In 

Clearly, (4.4) has the solution Qh = gh~^- Moreover, (4.5) is a system on mixed 
form admitting a unique solution provided that the finite element spaces satisfy the 

Babsuka Brczzi condition. However, the commuting diagram property satisfied 
by our finite element spaces immediately renders the Babuska-Brezzi condition 
satisfied (cf. Theorem A.4). □ 

5. Basic estimates 

In this section we establish a few estimates to be used later on, including squarc- 
intcgrability of the pressure and weak time-continuity of the density. However, we 
begin with the following lemma providing us with a renormalized formulation of 
the continuity scheme (3.1). 

Lemma 5.1 (Renormalized continuity scheme). Fix any m = 1,...,M and let 
(f?™''"™) ^ Qh >^ Vh satisfy the continuity scheme (3.1). Then (^/T''"™) 
satisfies the renormalized continuity scheme 

B{gt)4>h dx 

-AtJ2 f {B{8-){K ■ ^)+ + B{g^){u^ ■ ^)-) [</.,]r dx 

+ At [ b{g^)divu^<l>hdx+ [ B"{^{g^,g^-^)) [e^-'f <^h dx 
Jn Jn 



+ At^ [ B"($r(^-,^-))K]r('A^)-«-^)+ 

- B"{e{g-\ K]r (M+iK ■ dS{x) 

= j B{g^-^)(t>h dx, ^(t>h e Qh{^), 
Jn 



(5.1) 
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for any B € C[0, oo) n C^{0, oo) with B{0) = and b{g) := qB'{q) - B{g). Given 
two positive real numbers ai and a2, we denote by ^(01,02) and (,^{ai,a2) two 
numbers between oi and a-i (they will be precisely defined below). 

Proof. Since x ^ B' {g'^{x))(j)h{x) is piecewise constant, we can take B'{Q^)(j)f^ as 
a test function in the continuity scheme (3.3), yielding 

/ [e™-'] B\g^)^Hdx 
Jn 

= -AtJ2 [ {^+i< ■ + ■ '^r) B'{g^)<i>h dS{x) 

= -AtJ2 I {B'{g^)g^{K ■ y) - [Qh]aE • y)') <l>h dS{x). 

(5.2) 

A Taylor expansion yields 

B'{z){y -z) = B{y) - B{z) - B"{z*){y - z)\ 
for some number z* between z and y. Consequently, 

[QhU B'iQ"^) = [B{g^)]eE " B"{e''{e+. q'-)) [QhToe > 
W-'] B'ieT) = [B{g^-')]-B"{ag^,g^-')) [g^-']\ 

where 

e^^(e'|\ g"JKx) G [g'I^ix),g^ix)], x e dE, 

and 

e(^^^,e^')(^)eK"'(^).^>^^(^)], ^en. 

Inserting these identities in (5.2), recalling the definition of b, and applying 
Green's theorem, we achieve 



{BigT)-B{gl^-'))cl,^dx+ / B" i^eT , qT')) WrT 0h dx 
a Jn 

■■ -At [ b{g'^)diyu^(t)h dx 

- ^ / (B(e![')«" • + B{g^)iK ■ ^r) <t>h dS{x) 
EeEu JdE\dn 

+ AtJ2[ B'^e'^ig-^, g^)) [g^]l^ « • '^)-4>h dS{x). 



(5.3) 



EeEh 



Denote by / the second term on the right-hand side of the equality sign. Then, as 
in Remark 3.2, we have the identity 

I = AtJ2 [ {B{q+){K ■ ^)+ + B{g^){u^ ■ ^)-) dS{x). (5 4) 

rer^ •'^ 

Rec;alling that for each F G F/j we let and be the two elements sharing 
the face F and such that the normal component associated with F points from E+ 
to E- , we can write 



EeE^ JdE\dn 

= AtJ2 [ -B"{e^-{g^,g^))[gT]l{cl>H)-{<-yr ^^'^^ 
+ B"{e^^g^, g^)) [g^]l (</-,)+« • u)- dS{x). 
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Once we introduce into (5.5) the notations 

C (Q+,Q-)-=c. +(g+,g_), ^ {Q-,Q+)-=t. [0+,Q-), 
inserting (5.4), (5.5) into (5.3) yields the final result (5.1). □ 

In what follows we shall need a discrete Hodge decomposition. The following 
lemma is a consequence of (2.5). 

Lemma 5.2. Let {{gh,Wh,Ufi)}f^yQ be a sequence of numerical solutions con- 
structed according to (3.4) and Definition 3.1. For each fixed h > 0, there exist 
unique functions Ch e W°'-^ and e V^"'"^ such that 

< = curlCr + ^r, m = l,...,M. (5.6) 

Moreover, if we let Ch{t,x), Zh{t,x) denote the functions obtained by copending, as 
in (3.4), {Cr}m=i> {^r}m=i the whole o/(0,r] x SI, then 

Uh{t,-)^ curl Chi-, t) + Zh{-,t), te{0,T). 

We now state a basic stability estimate satisfied by any solution of the discrete 
problem given in Definition 3.1. 

Lemma 5.3. Let {{Qh,Wh,Uh)}f^yQ be a sequence of numerical solutions con- 
For any m = 1, . . . , M, we have 
P{Qh) dx 



structed according to (3.4) and Definition 3.1. For q > Q, set P{g) := t^Q^ ■ 



(5.7) 



m ,. 

+ E / P'\i{&l^eV)) {Qt't dx 

+ ^At / jw^f dx + ^At / jdivM^I^ dx 
fe=i "^^ fe=i "^^ 

+ ^^t \whf dx + ^At / jcurltu^l^ dx 

< / P(^o) dx + Cy^At / dx. 

Jq j.^-^ Jn 

Consequently, qh Eb i°°(0, T; L''(f])). 

Proof. Since P'{q)q — P{q) = p{s) and Qh > 0, it follows by taking = 1 in the 
renormalized scheme (5.1) that 

/ P(4) dx + At f p(^>^)div«^ dx+ [ P"{^{qIqI-')) [Ql-'f dx 
Ja Jn Jn 

+ Ai E / P"{e{e7^ e-)) [elt « ■ (5 8) 

- p'^eiQ"^, ^T)) [qI]Ik ■ ds{x) = [ p(/-i) dx. 

For k = 1, . . . , M and x G Urer^ 

k(.._i^^^{Q'lix),g'Lix)}, 1<7<2, 
'^^''>-\min{g1ix),g^_ix)}, 7 > 2, ^'-'^ 
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and note that 



rerl 

- P"{e{g'^, g+)) H] I « ■ dS{x), (5.10) 



>AtJ2 [ P"{e^)[et]l\ul,.\ dS{x). 



Next, by using = itjj as a test function in the first equation of (3.2) and then 
using the second equation of (3.2) with r]h = w^, we obtain the identity 

/ p(^^)divu^da;= (m+A) / |divu^|' dx+n [ \wl\^ dx- [ f'^u^dx. (5.11) 
JQ Jn Jn Jn 

Similarly, specifying Vh = curltu^ in the first equation of (3.2) yields 

/X / [curltuj^l^ dx = fll curl w'l dx. 
Jn Jn 

An application of Cauchy's inequality (with epsilon) then yields 

/ dx<C [ dx. (5.12) 

Jn Jn 

Thanks to (5.6), we can write ix^ = curl^^ + for with e W^'""" and 
^ '^h'^- Choosing rih = cml^h in the second equation of (3.2) gives 

/ |curlCfc|^ dx= w!^Ch dx<( {w^f dx) ( / |C^f dx) . 
Jn Jn \Jn / \Jn j 

Thus, since the discrete Poincare inequality (2.7) tells us that 

/ IC^I^ dx<C |curlCfc|^ dx, 
Jn Jn 

we arrive at the estimate 

/ |C^|% |curlCfc|^ dx<C [ \w^\'^ dx. (5.13) 
Jn Jn 

In view of the discrete Poincare inequality (2.6), we also have 



lz*=l 



dx < C I IdivixJ^I^ dx, 



which, together with (5.13), allow us to conclude 

/ \u^f dx= \Chf+\zh\^ dx<c( / \w^f dx+ |divu^|^ da;). (5.14) 
in Jn \Jn Jn J 

Now, by first inserting (5.11) into (5.8) and subsequently utilizing (5.10), (5.12), 
(5.14) and Cauchy's inequality (with epsilon) to treat the last integral appearing 
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/ 

JQ. 



in (5.11), wc acquire the estimate 
P{Ql)-P{Qt')dx 

+ / P"mlQ'k-'))[Qt't dx + MY, [ P"{g\)[e'l\l\ui-y\ dS{x) 

+ At dx + At / jdivti^l^ dx 

Jn Jn 



At 



I 

Jn 



Jn 



dx + At I \cm\wz\ dx<CAt 



I 

Jn 



dx. 



(5.15) 

□ 



Finally, summing (5.15) over k we conclude that (5.7) holds. 

The stability estimate only provides the bound j?(^?/i) L°°{0,T; L^{Q)). Hence, 
it is not clear that p{gh) converges weakly to an integrable function. Moreover, the 
subsequent analysis relics heavily on the pressure having higher intcgrability. In the 
ensuing lemma we establish that the pressure is in fact bounded in L^(0, T; L^(r2)), 
independently of h. 

To simplify notation, we denote the effective viscous flux by 

Pca{gh, Uh) = p{Qh) - (m + A) div Uh- (5.16) 
We will also continue to use this notation in the subsequent sections. 



Lemma 5.4 (Higher intcgrability on the pressure). LeA {(p/, , , tt;, )},j^q he a 
sequence of numerical solutions constructed according to (3.4) and Definition 3.1. 
Then 

p{Qh) &b Lmo,T)xn). 
Proof. For all m = 1, . . . , M, let e V^"'"^ be such that 

div< = Peff(^r,0 ^ 



Pes{gT,0 dx. 



Now, since the momentum scheme (3.2) gives 

Peff(^r> O div Vhdx = - I f^Vh dx, e V°-^, 

! Jn 



I 

Jn 



we can use as test function to obtain 



JjPM,K)f dx=^^ (^J^PM,0 dx 



fJTvi: dx. 



Hence, with e > 0, 



/ Peff(^r,0-^ / PeM^Odx 

= - / f;r< dx<^ [ \frf dx+e f \vi 

Jn Jn Jn 



(5.17) 



in 4e 
By the Poincare inequality (2.6) 



dx. 



I \v'^f dx<C [ 
Jn Jn 



Consequently, by fixing e small enough in (5.17), 



PM,K) dx 



dx<C 



dx. 



dx, 
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and thus 



^ \PM,0? ^ p {^j^PM,0 dx^ + C l/n' dx. (5.18) 
Now, due to the boundary conditions, 

/ Peff(er,0 dx= I p{q^) dx < C, 
Jn Jn 

where we also have put into use Lemma 5.3 and subsequently our assumptions on 

the source term / and the initial data qq. Hence, (5.18) allows us to conclude 



from which we obtain 



dx<C [1 



i/n dx 



m ■■ 



1,...,M, 



f; / \PM,Of dxKciT+f^^t I \fhf dx] 

m=l \ m=l J 

In view of the definition of Pcff , this immediately yields 

M „ 

J2 At / \p{Qh)f dx 

m=l ■'^ 

<c(T+J2^t [ |div<|' dx+Y^^tj l/ri' dx] 
which, due to Lemma 5.3 and / e i^((0,T) x f2), concludes the proof. 



□ 



We conclude this section by establishing a weak time continuity of the density 
approximation. For this purpose we shall need the following technical lemma, which 
provides a bound on the artificial diffusion introduced by the upwind discretization 
of the continuity equation. 

Lemma 5.5. Let {{gh,Wh,Uh)}i^yQ be a sequence of numerical solutions con- 
structed according to (3.4) and Definition 3.1. There exists a constant C > 0, 
depending only on the shape regularity of Eh, the size of \^\, and the final time T, 
such that 

i-T 

[Qh]aE i^h • i^)"(n^(/) - (j)) dS{x)dt 

ldE\dQ 



(5.19) 



< C'll^'/'ILoo(o,T;ioo(o))/l'. 

for any 4> e L°°{0,T; W'^'°°{n)). 

Proof. We shall need the auxiliary function 



B{z) 



7>2, 
7 < 2. 



Moreover, set 

</."(x) 



1 

At 



(s, x) ds, 



n' 



m 



1, 



,M. 



(5.20) 



Using B"{z) > for 2 > and Holder's inequality, we obtain 

M 

■ ~ ' K"WK"-^)"(C-'/'")rf5(x) 



E E At 

m=l EeEh 



dE\da 
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< 



X ( E E / (s"(er))"' K • H W - <i>"'\^ dS{x)\ 



=■■ h X h, 

where the "internicdiatc" numbers are defined in (5.9). 
If 1 < 7 < 2, then Lemma 5.3 can be apphed: 



Ii<C [ B{go) dx+y^At I 1/;^!' dx. (5.21) 
m=i 

However, (5.21) continues to hold in the case 7 > 2. This follows directly from the 
renormalized scheme (5.1), with (j)h := 1, together with the fact that 

M . 

V At / h{g'^)diYu'f^ dxdt 

^ (e ^^L i'^'^''"'' '^^j ' 



which is bounded by Lemmas 5.3 and 5.4. 
Next, using Lemma 2.9, we have that 



M 



E E At 



m=lEeEh •'dE\dn 
M 



{B"{g^))-'] dS{x) 



X E E At / K-H' dS{x) 

Vm=l EeE^ JdE\dQ ^ 



(5.22) 



Thanks to Lemma 2.11 and Lemma 2.10, 

[ \u^-i^f dS{x)<ch-^ [ lu'^f dx. 

JdE Je 

Moreover, since 

{B"{gJ^))-' < \g^ + g^f-'' < C{1 + g^ + g^), 
whenever 1 < 7 < 2, and ^i?"(^™)^ = ^, whenever 7 > 2, Lemma 2.11 



(5.23) 



also 



gives 



/ \{B"{gf)) dS{x)<Ch-^ [\E\+ f 

JdE '- \ Jeu 



X{E) 



leTr dx , 
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where Af{E) denotes the union of the neighboring elements of E. Observe that 



\ 



<h-^ 



Qh? dx 



\ 



U (i5UAr(E)) 



(5.24) 



where we have utihzed the fact that 



iN + 2) J^lg^f dx^ 



U (EUj^iE)) 



< (^ + 2)1^^1, 



EeEh 



IS 



which is true since the maximal cardinality of the set {{J\f{E) U E)r\ F} 
N + 2 hi any F e Eh- Inserting (5.23) and (5.24) into (5.22), we have arrived at 



h < /i' \\Dcj>\\ 



i~(0,T;i~(Q)) 

M 



h 2 



X T + 



^aW \eT\' dx] E^W K 

m=l J Vm=l ''^ 



dx 



< Ch ||-D^!>||j^oo(o,T;i=»(n)) ' 
where Lemmas 5.3 and 5.4 have been used to work out the last inequality. This 
concludes the proof of (5.19). 

To simplify the notation, let us introduce the interpolation operator 



□ 



{Ucf) it) = r-' + —^{r - r-'), t e (i™-\t™). 



(5.25) 



Lemma 5.6. Let {{Qh,Wh,Uh)}f^^Q be a sequence of numerical solutions con- 
structed according to (3.4) and Definition 3.1. Then 

d 



dt 



{UcQh) L^(0,T; W-^'\^)). 



Proof. Fix (t) e L°°(0, T; W^^°°{n)), and recall the definitions of (/>™, (j)"}^, cf. (5.20). 
The continuity scheme (3.1) with (/)™ as test function reads 

{^cQh)cr dxdt 

(5.26) 



At 



dt 



= AtY, [ ■ + Q7i< ■ ^r) [fflr dS{x). 

Since the traces of taken from either side of a face are equal, we can write 
E / {e-{<-yr + Q':^{u^-u)-)m^dx 

= E / • + K" • ^)+) [c - n ds{x), 



E 



EeE^ JdE\dn 



{g^iu^ ■ u)+ + Q^{u^ ■ v)-) (C - r) dS{x), 
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Jn JdE\dn 
= [ eh< ■ Dr dx+Y. I ^Sh]oE « • - r) dS{x). 

JQ. EeEh 

By summing (5.26) over m, taking absolute values, and using the above identity, 
we find 

M , 

J2^At j^-{ncQh)<r dxdt 

M „ 

At / g^KDcP"^ dx 

M . 

+ YT. [g^],^ « • '^)-(C - r) dS{x) 

rn=l EeEh J dE\dn 

Using Lemma 5.5, together with an application of Holder's inequality, we deduce 

M 



m 1 O 



M \ 2 / M \ 2 

Y At / \qT\' dx) E A* / \<\' \\mL^iO,T;L^m) 
m=l •'^ I \m=l •'^ J 



< 



T;I,~(n)) • 



By Lemmas 5.3 and 5.4, the first two factors on the right-hand side is bounded, so 

we conclude that 



/ / ^ i^cQh) <^ dxdt 
J At Jn "t 

M . , 

VAt / - (ncQh)'rdx 



<C(l + /i5)||Z?</,||^„(o_^.^„ 



□ 



6. Convergence 

Let {{Qh^Wh,Uh)}f^yQ be a sequence of numerical solutions constructed accord- 
ing to (3.4) and Definition 3.1. In this section we establish that a subsequence 
of {{QhiWhTUh)}h>o converges to a weak solution of the semi-stationary Stokes 
system, thereby proving Theorem 3.3. The proof is divided into several steps: 

(1) Convergence of the continuity scheme. 

(2) Weak sequential continuity of the discrete viscous flux. 

(3) Strong convergence of the density. 

(4) Convergence of the velocity scheme. 

Our starting point is that the results of Section 5 assure us that the approximate 
solutions {wh,Uh, Qh) satisfy the following /i-independent bounds: 

Qh Gb L°°(0,T;LT(n))nL27((o,T) xfi) 



and 



div,2/ 
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Consequently, wc may assume that there exist functions g,w,u such that 
f?/, mL°-{0,T;L^{n))nL^^{{0,T)xn), 
w^'-^^w, inL^{0,T;Wr'-'m), (6-1) 
Uh'^^^u, in L\0,T;Wt'\^))- 

Moreover, 

where each ''^'^ signifies weak convergence in a suitable space with p > 1. 

Finally, Qh, Qh log Qh converge respectively to g, glog g in C([0, -^weak( ^)) fo'' 
some 1 < p < 7, cf. Lemma 2.2 and also [9, 16]. In particular, g, glogg, and ^log^ 

belong to C{[0,T];Ll^^^{n)). 

6.1. Density scheme. 

Lemma 6.1 (Convergence of gnUh)- Given (6.1), 

QhUh Qu in the sense of distributions on {0,T) x Ct. 
Proof. By virtue of Lemma 5.2, there exist sequences {Ch}h>Oj {zh}h>o satisfying 

Uh{-,t) = curlCft(-,t) + Zk{-.t), 

Ck{;t)eW'^^^, ZH{t,-)€V^'^, 

for all t e (0,T). In Lemma 6.2 below we prove that 

curio,, ^ curie in L^{0,T; L^{n)). 

As a consequence, curl (^h Qh curl C, g the sense of distributions. 
It remains to prove that 

QhZh QZ in the sense of distributions. 

To this end, we adapt the proof of [16, Lemma 5.1] to our specific discrete setting. 
We begin by introducing the regularized field = ★ z/j, where /t^ is a standard 

(a:) 

regularizing kernel and -k denotes the convolution product (in x). Lemma 6.3 

{x) 

guarantees that 

114 - ^ft|li,2(o,T;i2(n)) ^0 as e ^ 0, uniformly in h. 

In addition, for any k and p, since z^ G i^(0, T; W^-''p{U)) we have that z^ z^ 
in L2(o,T; W*='P(n)). Moreover, z<= z in L2(0, T; ^^(f^)). Hence, by writing 

QhZh = ^/i(^/i — + ^/j^;^ it sufiices to prove ghZ% ^° £»z^ for each fixed e > 0. 
Next, let us introduce auxiliary functions Z^'^ G Vh, m= 1, . . . , M, defined by 

m 
k=0 

We extend {^^"1™=! a function defined on {-At,T] xCthy setting 

z^(t,.) = zr(-), te(t'"-sn. m=i,...,M, 

and Z^(t, •) = Zf, for t G {-At, 0]. In view of the regularity of z^, 

Z^(t, •) ^ Z'=(i, •) = / z'(s, •) ds in C'=(n) for any A; > 0, (6.2) 
Jo 

uniformly in t on [0,T]. 



24 



K. H. KARLSEN AND T. K. KARPER 



Now, we write 

'7t,m m-l i^e,m-l m 

Qh^h - ' 
which alternatively can be written as 

d d 

on tm] X 17, m = 1, . . . , M. 

Fix (j) G C^((0,r) X n). Summation by parts gives 

/ ^ncighZf,)<t>dxdt 
Ja 'J^ 

= - I I Bhit - At,x)Zl{t ~ At,x)- aicM dxdt, 
J At Jn 

where (phit, •) = St It^--^ '?^(^' ") * ^ (im-i,tm)- 

Thanks to (6.1) and (6.2), QhZf, qZ" in L'^^{0,T; L'^^{n))nL°°{0,T; ^^(fi)), 
and hence 

(ghZ^) — ^ — {qZ'^) in the sense of distributions on (0, T) x fl. 
In addition. Lemma 5.6 tells us that ^ O^cQh) Gfc W^^'^{^1)), and thus 

z^(.-A,.)|(n.,.)''^°z^|, 

in the sense of distributions on (0, T) x Cl. 

We conclude observing that gz" = ^ {gZ^) - Z"^. □ 

In the proof of the previous lemma we utilized 

Lemma 6.2. Given (6.1), define {(Cfn ^/i)}/j>o terms of the decomposition 

Uh{t,-) = cmlCh{t,-) + Zh{t,-) with Chit,-) e W°'^, Zh{t,-) e V^'^, t e (0,T). 
Then 

Wh w, curia curie L^{0,T; L^{n)). (6.3) 

Proof. Subtract the first equation of (3.2) with Vh = curl^™ from fi times the 
second equation of (3.2). Multiplying the result with At and summing over all 
m = 1, . . . , M yields 

/ fi curl rill curl C/» ~ A* curl curl dxdi; 
Jo Jn ^g_4^ 

/ ^J'Whr]h - fh cml^h dxdt, 

Jn 

for all r]h,^h that are piecewise constant in time with values in Wft,(r2). Fixing 
r/,^ e C^((0,T) X f2), we use in (6.4) the test functions 

^ft(i,-) = C(-)-^_^^ ^n^^(.,s) ds, ie(t™_i,U,'«=i,---,^- 

a(i,-) = C(-):=^^^ ^nfC(-,s) ie(i™_i,i„), m = l,...,M. 
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// 

Jo Jn 



Due to Lemma 2.9, curl^?i curl^ and c\xv\r]h cmlrj in L^(0, T; L^(r2)). As a 
consequence, keeping in mind (6.1), we let /i — > in (6.4) to obtain 

rT p 

jjL curl r] curl C — A* curl w curl ^ dxdt 

T (6-5) 
= / I fiwT)- fcmlC dxdt, Vtj,C e C~((0,T) X Q). 

Since C~((0,r) x Q) is dense in ^^(o^ T; W"o"''^(Sl)) ([12]), wc conclude that 
(6.5) holds for all 77,^ e L'^iO,T; Wo"''^(f^)). Hence, taking 77 = to, ^ = C in (6.5), 

0= / / n\wf - fcml^ dxdt. (6.6) 
Jo Jn 

Next, setting r]h = "Wh and = in (6-4), we observe that 

r-T 



= / jj, \whf - fhcuil^h dxdt. 
Jo Jn 



10 Jn 

Letting /i — > and comparing the result with (6.6) reveals that 

lim / / n\whf' dxdt = / / dxdt, 
'^^ojo Jn Jo Jn 

which implies the first part of (6.3): 

Wh^w in L^{0,T;L^{n)). (6.7) 

To prove the second part of (6.3), we make use of rjh = as a test function in 
the second equation of (3.2), sum the result over m = 1, . . . , M, and subsequently 
send h to zero: 

lim / / IcurlCfel^ dxdt = lim / / WhCh dxdt 
h^oJo Jn '^^Wo Jn 

'^==^ / / 'wC dxdt = I I I curl^l^ (ia;rff, 
Jo Jn Jo Jn 

where the last equality follows by arguing along the lines leading up to (6.6). □ 

During the proof of Lemma 6.1 we made use of a spatial compactness property 
stated in the next lemma. 

Lemma 6.3. Given (6.1), define {{Ch,Zh)}f^^Q in terms of the decomposition 

uh{-,t) = cMTKH(t,-) + zh{,-,t) With a(-,t) e w°'^, zh{t,-) e V°'^, forte (o,r). 

Then, for any ^ € M.^ , 

\\zhit, •) - Zhit, ■ - mlHo,T;mn,) < C (iCl"^ + ||div^^||i.(o_y.^.(j,)) , 
where fi^ = {a; € f2 : dist(a;, dO) > 

Proof For each t € (0,T) we know that Zfc(i, •) e Vj[''"'"(f2), so Theorem A.l can 
be apphed to give 

\\zn{t, •) - Zf,{t, ■ - Ofmn,) < C (l^l"^ + 1^1^) lldiv Zft||i.(^) , 

where C > is independent of h,^,t. We conclude by integrating over (0,T). □ 

Lemma 6.4 (Continuity equation). The limit pair {g,u) constructed in (6.1) is a 
weak solution of the continuity equation (1.1) in the sense of Definition 2.5. 
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Proof. Fix a test function 4> G C^ ([0, T) x fi), and introduce the piecewise constant 
approximations 4>h n^t/i, c/i™ := Tl'^4>™, and ^™ := /jL-i •) rft- 

Let us employ (/>™ as test function in the continuity scheme (3.1) and sum over 
m = 1, . . . , M. The resulting equation reads 

M 

TO=1 ^ 

M 

= J2T.^t {g^{uT-u)+ + g^{uT.u)-)[<j>^]^ dS{x). 
As in the proof of Lemma 5.6 we can rewrite this as 

M 

= ^ At / p;r<^'A™ ^^2; 

m=l •^f' 

M „ 
E'^ir ™_i JdE\dn 



(6.8) 



BeBfc m=l 

/ QhUhDcj) dxdt 

T 



+ E / / • {(t>h- (!>) dS{x)dt. 

w.^w.. Jo JdE\dn 



~(0,T;I,~(n)) • 



Lemma 5.5 tells us that 

E / / MaB(M?i-I^)"((/>ft-(/>) ^^(x)^^ 

In view of Lemma 6.1, 

lim / / ghUh,D(j) dxdt = / QuD(j) dxdt. 
f^^^Jo Jn Jo Jn 

Summation by parts gives 

^ At / - {UcQh) C rf^rf^ 

TO=1 ^ 

= - / / Qh{t-At,x)^ (Uc(l>h) dxdt - [ eUidx 
J At Jn ot Jq 

— / gcpt dxdt — / go<j>{0, x) dx. 
Jo Jn Jn 



where (6.1), together with the strong convergence g^ go, was used to pass to 
the limit. Summarizing, letting — > in (6.8) delivers the desired result (2.2) □ 

6.2. Strong convergence of density approximations. The instrument used 

to establish the strong convergence of the density approximations is a weak 
continuity property of the quantity Pes{gh,Uh) defined in (5.16). To derive this 
property we exploit our choice of numerical method and the boundary conditions; 
specifically, the finite element spaces, which are chosen such that (6.10) below holds. 
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Lemma 6.5 (Discrete effective viscous flux). Given the convergences in (6.1), 

lim / / Peff{gh,Uh) Qh dxds = / / Peff{Q,u) g dxds, Vte(0,T). 
Proof. For each m = 1, . . . , M, consider the problem 

div = -^Jjh dx, C ^= " ^ £ _^ dt, (6.9) 

where gJJ = — go- Observe that JqQ^^ dx = 0. Indeed, using the continuity 
scheme (3.1) and the continuity equation satisfied by the limit g, cf. Lemma 6.4, 

/ gh dx = / g1 dx, / H^^) dx = / g dx = / go dx 

JQ Jfi Jn Jn Jo. 



Thus, there exists a unique sohition G V^^'^ of (6.9). We denote by Vh{t,-), 

qh{t,-) the usual "piecewise constant" extensions of {w/^'j^^i, {(/"j^^i *° (0'^)- 
Utilizing as test function, the velocity scheme (3.2) reads 

/ PM,K)C dx=^ ( qldx(( PM,K) dx) - [ /r^r dx. 
Jq I"I Jn \Jq / Jn 

(6.10) 

Multiplying by At, summing over m, and using the definition of q^, we arrive at 

/ / Pes{Qh,Uh){0h - e) dxds = [ Qhdxl [ [ Pe«{gh,Uh) dxdt) 
Jo Jn l"l Jn \Jo Jn J 

\ fhVh dxds, 
Jn 

for any t e (0,T). 

In view of Theorem A.l, we have that Vh ^ in L^(0, T; L^(fi)). Since fh^f 
in L^((0, T) X ri) and gJJ dx 0, we conclude the desired result 



lim / / Pef{{gh,Uh){gh- q) dxds = 0. 

h^OJo Jq 



□ 



We are now in a position to infer the sought-after strong convergence of the 
density approximations. 

Lemma 6.6 (Strong convergence of gh)- Suppose that (6.1) holds. Then, passing 
to a subsequence if necessary, 

Qh ^ Q a.e. in (0,T) x fl. 

Proof. In view of Lemma 6.4, the limit {g, u) is a weak solution of the continuity 

equation and hence, by Lemma 2.7, also a renormalized solution. In particular, 

(g log g)^ + div ((o log g) u) = g div u in the weak sense on [0, T) x 17. 

Since t g^ogg is continuous with values in some Lebesgue space equipped 
with the weak topology, we can use this equation to obtain for any t> Q 

{glog g) (t) dx — / gologgodx = — / / gdivudxds (6.11) 
Jn Jo Jn 

Next, we specify t/i^ = 1 as test function in the renormalized scheme (5.1), 

multiply by At, and sum the result over m. Making use of the convexity of z log z, 

we infer for any m = 1, . . . , M 

/ g^ log g]^ dx- [ gl log gl dx < -J^At I g^ div dxdt. (6.12) 



Jn 



f glloggl dx<-y^At j g'f:dWu^ dxdt. 
Jn Jn 
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In view of the convergences stated at the beginning of this section and strong 
convergence of the initial data, we can send /i —> in (6.12) to obtain 

/ [B^og Qj{t) dx — / Qologgodx<— / / gdivudxds. (6.13) 
Jn^ ' Jn Jo Jn 

Subtracting (6.11) from (6.13) gives 

J (^glog Q — glog gj (t) dx < — J J gdivu—gdiyudxds, 
for any t G (0,T). Lemma 6.5 tells us that 



gdivu— ^divu dxds = / / — g'^g dxds > 0, 

Jn M + JO in 

where the last inequality follows as in [9, 16], so the following relation holds: 



glogg= glogg a.e. in (0,T) x ft. 
Now an application of Lemma 2.1 brings the proof to an end. □ 
6.3. Velocity scheme. 

Lemma 6.7 (Velocity equation). The limit triple {w,u,g) constructed in (6.1) is 
a weak solution of the velocity equation (1.2) in the sense of (2.3). 

Proof. Fix {v,r]) G C^((0,T) x Q), and introduce the projections Vh = H^v, 

Wi = r? and = it It^-i dt, = Vh d.t. 

Utilizing v"^ and 77™ as test functions in the velocity scheme (3.2), multiplying 
by At, and summing the result over m, we gather 

/ n curl WhVh + [(/X + A) div Uh - p{gh)] div Vh dxdt = / fhVh dxdt, 
Jo Jq Jo Jq 

WhVh - Uh cml rth dxdt = 0. 

(6.14) 

From Lemma 2.9 we have Vh >° v in L^(0, T; Wq '^'^(Q)) and ijh ^" V in 
L^(0, T; W,J"'''^(ri)). Furthermore, by Lemma 6.6 and the first part of (6.1), 

p{6h) p{q) in -^^((0,7") x Ct). Hence, we can send /i — > in (6.14) to ob- 
tain that the limit {w,u, g) constructed in (6.1) satisfies (2.3) for all test functions 
{v, T]) e C^{{0, T)xn). Since C^HO, T) x Q) is dense in both ^^(o, T; Wo^'^{^)) 
and L2(o, T; Wo"''''^(f})) [12], this concludes the proof. □ 



Appendix A. Compactness of functions in V^'^{Vl) 

In this appendix wc prove that discrete weakly curl free approximations in 
V°'-'-(r2) with L"^ bounded divergence possesses an L"^ space translation estimate, 
which was previously needed to conclude the weak convergence of the product 
ghUh to the product of the corresponding weak limits gu. As part of the proof, in 
Lemma A. 5 we show that if a sequence {vh}h>o belongs to V^"^{Vt) and besides 
satisfies divvh &b L'^i^), then {curlt;/i}?i>o and {diyvh}h>o are actually compact 
in W~^''^{il). Thus, strong _L^(f2) convergence of a subsequence of {fh}/i>o follows 
directly from the div-curl lemma. However, this is not sufficient to conclude the 
sought after convergence of gnUh (cf. Subsection 6.1). The problem is a lack of 
temporal control of the velocity approximations Uh- 
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A.l. Space translation estimate. The argument is inspired by Brenner's work 
[2] on Poincare-Friedrich inequalities for piecewise vector fields. The basic 
idea is to project the relevant function into the Crouzeix-Raviart element space 
and then use the standard translation estimate satisfied by functions in this space 
(cf. Stummel [21]). Then, since the relevant function is discrete weakly curl free, 
we have sufficient control on the curl to suitably bound the projection error. 

The Crouzeix-Raviart element space is defined as a non-conforming Pi element 
for each component of the vector field. That is, on each element E g the 
Crouzeix-Raviart polynomial space R{E) is given by 

R{E) = [¥\E)]^, 

where Pi (E) is the space of linear scalar fields on E and N is the spatial dimension. 

The degrees of freedom of R{E) are the average integrals over the faces of E. The 
Crouzeix-Raviart element space, denoted Rh{^), is formed on Eh by matching 
the degrees of freedom on each face F e F/j. Hence, Rhi^) is discontinuous across 
element faces and thus leads to non-conforming discretizations of H^{il). However, 
it has the property that for any Vh € Rh, /r[^/i] dS{x) = for all F e Th- 

Theorem A.l. Given Zh G V^^''^! there exists a constant C > 0, depending only 
on and the shape regularity of Eh, such that for every vector £, G M.^ , 

\\zhi-) - Zh{- - OWmn,) < C (l^l"^ + 1^1') ' ||divz,,||^.(^) , (A.l) 
where f2j = {a; € : dist(a;,9ri) > ^}. 

Proof. Let us introduce an interpolation operator ; Vh —> Rh by specifying 
n^Zh dS{x) = 4 / i^h} dS{x), VF e F^, 



1 



where {•} denotes the average of the traces from the two sides of F. According to 
Brenner [2], we have the following error estimate: 



\Zh ■ 



By a standard decomposition of vector fields, 

[zhlr = i(zh ■ t^)Hr ~ iizh X 1^) X - [{zh x i^) x v] 



(A.2) 



r ' 



where the last equality follows for the reason that [zh ■ i^]r = for all F e F^. To 
have the above decomposition well-defined in two dimensions, we set {zh xy)xu = 

— {zh X v)t, where t is the tangential vector. 

Since v is constant on each F e F?i with \v\ — 1, 



[zh] dS{x) 



J [zh X I/] dS{x) 





2 




X v 


< 











[zh X u] dS{x) 



VF e F^. 



As a result, applying (A. 11) of Lemma A.2 below to (A.2) yields the error estimate 



\zh-n^Zh\\l,^^^<clYl h'WDzhWl.^^^ 
\EeEh 



4-N 

- h 2 



jdiv Zh\ 



L2(0) 



(A.3) 
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To continue, fix an arbitrary ^ e . By the triangle inequality, we write 

<c(||z,(.)-n«z,(.)||l.(^^) 

+ ||n«z,(.)-n^z,(.-o"' 



(A.4) 



+ \\U^ZH{--0-Zh{--0\ 



which transfers the translation onto the projected function H^z. 

Since II^z/j is a function in the Crouzeix-Raviart element space, Stummel's work 
[21, Theorem 2.1] can be applied, yielding 



|n^z,(-)-n^z„(--0| 



<c{h' + \e) J2 ll^n«z,||^,(^) 

EeEh 



(A.5) 



where the constant C only depends on f2 and the shape regularity of E^. 
Utilizing (A.3) and (A.5) in (A.4) gives 



<C{h' + \^\') E \\Dz,\ 
EeEh 



• h^" ||divz,,||^2( 



Since \Dzh\ = |div Zh\ on each E & Eh, this immediately yields 

\\zh{-) - zh{- - m^n,) < C {h^^ +h'+ ICI') ' l|divz;,||^.(^) 



(A.6) 



Finally, let us argue that (A.6) implies (A.l). To this end, fix any ^ € and 
h > 0. There exists a shape regular partition Gh of Q into triangles/tetrahedrals 
such that LiEeGh^ — ^EeEu^i ^ach E G Gh has a non-empty intersection with at 
most one element E G Eh, and such that 

max diam(i!;) < . 

E£Gh 3 

Next, let V|^|(0) denote the first order div conforming Ncdclcc clement space of 
first kind formed on the mesh G^,, and let n|^| : Vh{^l) V|j|(ri) denote the usual 
projection into this space. 

Now, for any Zh S V^''^{fl), we calculate 



\\Zh{-) - Zhi- - OWl^Q^) 



< 



z,(-)-n[;|Z^(-) 



+ 



uLzhi- - - - 



nf^|Z^(-)-nf^|Z^(--0 

2 



2 



(A.7) 



By Theorem A.l, 

KM-)-'nyh{--0 



< 



C(kl'^ + kl')'||divz,||^.(^). (A.^ 
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By Lemma 2.9, 

z^(-)-nf^|Z^(-) 



E&Eh 



(A.9) 

□ 



Inserting (A. 8) and (A.9) into (A. 7) completes the proof of (A.l). 

A. 2. Tangential jumps. In the proof of Theorem A.l we harnessed 

Lemma A. 2. Given Zh € V^'^ , there exists a constant C > 0, independent of h, 
such that 



j [zh X v] dS{x) 



<Ch^ ||div;s^||^. 



(n) 



vr e r 



h, 



and 



S}1 



Zh X v] dS{x) 



<C\n\^ ||divzfc||^2 



(Q) 



(A.10) 
(A.11) 



Proof. Let <p G Wg^'^(n). In virtue of Lemma A. 5 below, 



Zh curl (f) dx 



<Ch\\cl>\\ lldiv ZhWh'^ 



2(a) 



Applying integration by parts (2.4), keeping in mind that c\ix\zh\E = for all 
E e Eh, yields 

2_] / 4>[zh X v\ dS{x) = / Zhcml<p dx, 



and so 



V [ (P[zh X u] dS{x) 



<chu\\ lldiv 

■^/l II 1,2 



2(a) 



(A.12) 



The bound (A.12) serves as the starting point for proving (A. 10) and (A.11); 
the remaining objective is to construct a suitable test function 0. Fix T € F^. Let 
E-, E^ denote the two elements in Eh sharing the cgde/face T, where E-, E+ are 
chosen so that points from E- to E+. In view of Lemma A. 3, we can choose a 
continuous piecewise linear (scalar) function ^ on F such that 



Ijf dx=j^l^fdx. 



v/€Pi(r). 



Denote by cjjdE the extension by zero of cp to (9_E_ \ F) (J (9S+ \ F), and fix a 
piecewise affine function (pE on E- U E^ such that 4'E\g^_^^JQ^ = 4'dE- Clearly, 
(pE can be chosen such that 

I I < C/i"^ in the interior of E_UE+. 

Finally, let (f)^ denote the extension by zero of i^^; to all of Q.. 

The function (^r possesses the following properties: (j!)r € Wo'^(fi), (;&r|f = for 
all f e F,, such that f F, and 



11-^: 



H'i.2(£;_) 



wi.2(£;+) 



< C/i^ (1 + h-^) . 



li N = 2 (curl is scalar), then we opt for = (/)r in (A.12) to obtain 



L 



Zh X v] dS{x) 



N 

< Ch^ \\divZh\\L2^a) 



(A.13) 
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If A'' = 3, (A. 13) still holds. Indeed, to conclude we can in (A. 12) successively take 
= [0p,O,O]'^, (p = [0, (^r,0]^, and = [0,0,(j)r]'^ . Since T e T/j was arbitrary, 
this concludes the proof of (A. 10). 

To establish (A. 11), we introduce the test function 

(Pr = <PrJ^[zhXiy]dS{x), VF e T^, 

where is constructed as above with the additional requirement that 
supp 4>r P) supp (/>f = 0, Vf 7^ r. 

We have 

sup|£)0r| < Ch-^ 
and, for each F e Th, 

^ (f>rf dS{x) = ^ i^h X H dSix)^ (^^ / dSix)j , V/ G ¥\T). 
Finally, we set <p := X^rer^ '/*r; this function satisfies (p G Wo^'^(f2) and 

sup < C/i^^ max 
xen rer^ 

The last inequality follows from (A. 10). A direct calculation gives 

||£'0IL2(n) < Ch^ ||divz,,||^2(a) • 
Setting <p as test function in (A. 12) immediately gives the estimate 



which is (A. 11) 




The next lemma provides us with the specific test function that was brought into 
service in the above proof. 

Lemma A. 3. Fix any F G F^. There exists a continuous piecewise linear (scalar) 
function (j) onV such that (j)\dT = 0, \4>{x)\ < 1 Va: G F, and 

j^4>f dx=^ jj dx, V/GPi(F), (A.14) 

where N is the spatial dimension. 

Proof. Let b denote the barycentric middle point with respect to the vertices of 

F. Let T/i be the triangulation of F obtained by setting 6 as a vertex in addition 
to the vertices of F. On let LhiT) denote the standard finite element space 
of continuous piecewise linear functions. Any function (j)h G LhiV) is uniquely 
determined by it's value at the vertices. 

Now the relevant test function 4> G ^/i(r) is obtained by requiring 

4){b) = 1 and (t>{vi) — 0, for all vertices Vi at dO.. 

By direct calculation it can be verified that p satisfies (A.14). □ 



L 



Zh X v] dS{x) 



VF G F 



ft' 



[zh X v] dS{x) 



<Ch^ lldivz^l 



L2(n) 
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A. 3. Negative space compactness of the curl. In the proof of Lemma A. 2, 
the essential ingredient was an estimate on the W~^'^ norm of cml Zh- In this 
subsection, we prove this result. 



(A.15) 



Theorem A. 4. Consider the mixed Laplace-type problem 

curl w — D div u = f, w = curl u in fl, 
u-i' = 0,wxi' = on dCl, 

where we assume f G L'^{Cl). There exists a pair 

satisfying (A.15) in the weak sense. Moreover, there exists a pair 

{wh,Uh)eWh{n)xVh{n), 

satisfying the corresponding mixed finite element formulation of (A.15). Finally, 
the following error estimate holds: 

W-w - ^hh^n) + 11-" - UhWv^ < Ch'WfhHQ), (A.16) 

where the convergence rate s G [1/2,1) depends on the regularity of dQ. If dQ is 
Lipschitz and convex, (A.16) holds with s = 1. 

Proof. For example, cf Theorem 7.9 in [1]. □ 

Lemma A. 5. Let {zh}h>Q be a sequence in V^''^ for which ||divz/(||j;_2(Q) < C, 
where the constant C > is independent of h. Then 

||curlzfc||^_i,2(n) < C/i ||divzfc||^2(a) , 

for some constant C independent of h. 

Proof. To prove this lemma, we will use the mixed system (A.15) to define a new 
operator. To motivate the construction, consider the problem 

- A0 = curl in n, • z/ = 0, curl d xu = on dCl, (A.17) 

for some given e Wq^'^^'^ {^l) . By utilizing DA~^ div 9 as test function in the 
weak formulation of (A.17), where is the Neuman Laplace inverse, it is easily 
seen that the weak solution 9 of the system (A.17) is divergence free. Furthermore, 
we can set w = cml 9 and integrate by parts to conclude that the pair {w,9) is 
also the unique weak solution of the mixed Laplace system (A.15) with / = curl 0. 

Now wc define a new operator II'* : Wq^"^'"^ Wh as the unique function 
Il^(f) e W/j satisfying the finite element formulation: 

/ cmlCn^ (p)vh + div 9hdivvh dx = / cuvl(pvh dx, Vvh^Vh, 
Jn Jn ^^-^g^ 

h 

a 



curlr]h9h dx ^ I {U"(f))r]h, dx, Wr]h e Wh- 



The existence of such a function H'^cf) is given by Theorem A. 4. Using the fact that 
div 9 — 0, the error estimate (A.16) yields 

||nV - curl0||^,(^^ + ||0^ - 0||i2(n) + ||div0,||^2(a) < Ch^ I|curl0||^2(o) • (A.19) 
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Let Zh be as stated in the lemma. Since Zh is orthogonal to functions in Wh-, 

||curlz,,||^_i,2(n) = sup J— 

Zh curl(0 — n''0) dx\ 



= ii^ii 

(A. 18) I Jo div Zh div 6h dx I 

= sup 



0ewj'^ ll'/>llwi.2(Q) 
< C/i ||divzfc||^2(n) , 

where we have used (A. 19) to derive the last inequality, specifically the estimate 

I|div0/j||j,2(n) < C/i||curl0||^2(n) < C'/i ll0llvKi.2(n) • 
This concludes the proof. □ 
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